function M = madist(mdl, X, op)
%MADIST Computes the Mahalanobis distances to the Gaussian center
%
% [ Syntax ]
%   - M = madist(mdl, X)
%   - M = madist(mdl, X, 'shared')
%
% [ Arguments ]
%   - mdl:      the model or an array of models
%   - X:        the sample matrix
%   - M:        the resulting matrix of distances
%
% [ Description ]
%   - M = madist(mdl, X) computes the distances from the samples in X
%     to the center of Gaussian model(s).
%
%     Suppose there are n samples, then X should be a matrix of size
%     d x n, with each column giving a sample. 
%
%     If mdl is a single model object, then in the output, M is a row
%     vector of size 1 x n. If mdl is an array of K models, then M is
%     a row vector of size K x n. 
%
%   - M = madist(mdl, X, 'shared'). If you know that all models in mdl
%     share the same covariance matrix, you can use this syntax to tell
%     the function, such that it can make use of this fact to accelerate
%     the computation.
%
% [ Remarks ]
%   - If mdl is an array of multiple models, then all models should be 
%     in the same cache level. 
%
% [ History ]
%   - Created by Dahua Lin, on Dec 24, 2007
%

%% parse and verify input

assert(isnumeric(X) && ndims(X) == 2, ...
    'sltoolbox:slcdiaggauss:madist:invalidarg', ...
    'X should be a numeric matrix');

K = numel(mdl);
d = getdim(mdl(1));

assert(size(X,1) == d, ...
    'sltoolbox:slcdiaggauss:madist:invalidarg', ...
    'The dimension of samples does not match that of the models.');

if K > 1
    assert(isequal(mdl.clevel), ...
        'sltoolbox:slcdiaggauss:madist:invalidarg', ...
        'all models in mdl should be in the same cache level');
end

clevel = mdl(1).clevel;

if nargin >= 3
    assert(ischar(op) && strcmp(op, 'shared'), ...
        'sltoolbox:slcdiaggauss:madist:invalidarg', ...
        'the 3rd argument should be a string "shared"');
    
    shared = true;
else
    shared = false;
end


%% skeleton

switch clevel
    case 2
        M = compute_cl2(mdl, X, shared);
    case 1
        M = compute_cl1(mdl, X, shared);
    case 0
        M = compute_cl0(mdl, X, shared);
end


%% core computation routines

% Level 2

function M = compute_cl2(mdl, X, shared)

if shared || numel(mdl) == 1
    invsm = mdl(1).invsm;
else
    invsm = [mdl.invsm];
end
    
M = compute_core(X, numel(mdl), shared, ...
    invsm, [mdl.ism_mu], [mdl.mu_ism_mu]);

% Level 1

function M = compute_cl1(mdl, X, shared)

if shared || numel(mdl) == 1
    invsm = 1 ./ mdl(1).sigma;
    ism_mu = bsxfun(@times, invsm, [mdl.mu]);
else
    invsm = 1 ./ [mdl.sigma];
    ism_mu = invsm .* [mdl.mu];
end

M = compute_core(X, numel(mdl), shared, ...
    invsm, ism_mu, [mdl.mu_ism_mu]);

% Level 0

function M = compute_cl0(mdl, X, shared)

if shared || numel(mdl) == 1
    invsm = 1 ./ mdl(1).sigma;
    ism_mu = bsxfun(@times, invsm, [mdl.mu]);
else
    invsm = 1 ./ [mdl.sigma];
    ism_mu = invsm .* [mdl.mu];
end

mu_ism_mu = sum([mdl.mu] .* ism_mu, 1);

M = compute_core(X, numel(mdl), shared, ...
    invsm, ism_mu, mu_ism_mu);


% common computational skeleton

function M = compute_core(X, K, shared, invsm, ism_mu, mu_ism_mu)
% inv_sm: d x K
% ism_mu: d x K
% mu_ism_mu: 1 x K

% ord-0 term
t0 = mu_ism_mu';

% ord-1 term
t1 = ((-2) * ism_mu)' * X;

% ord-2 term
t2 = invsm' * (X .* X);

% combine
M = combine_terms(t0, t1, t2, K, shared);



%% supporting functions

function M = combine_terms(t0, t1, t2, K, shared)

if K == 1
    M = t1 + t2 + t0;
    
elseif shared
    M = bsxfun(@plus, t1, t2);
    M = bsxfun(@plus, M, t0);
    
else
    M = t1 + t2;
    M = bsxfun(@plus, M, t0);
    
end

